version 16
set more off
capture log close
cls
clear

*	Set color scheme
set scheme s2color

* 	Use the macro $S_DATE
global date : di %tdCCYY.NN.DD date("$S_DATE","DMY")
cd "/Users/davidflood/Box Sync/Diabetes systematic review"
log using "Log files/log_${date}.log", replace

import excel "Extraction/2020_10_26 - Extraction.xlsx", sheet("Sheet1") cellrange(A1:BK41) firstrow

**********************************************************************
*																	 *
*   SYSTEMATIC REVIEW AND META ANALYSIS OF DIABETES HEALTH SYSTEMS	 *
* 	INTERVENTIONS													 *
*																	 *
**********************************************************************

/*******************************************************************************
Part 1: Set up dataset
*******************************************************************************/

//	Drop irrelevant variables

drop criteria total_participants_string comments author_emails acknowledgements

//	Clarifying imported variables

*	Descriptive variables

label variable author_year "Author and year"
label variable study_title "Study"
label variable study_design "Study design"
label variable country "Country"
label variable wb_income_group "World Bank income group"
label variable wb_region "World Bank region"
label variable epoc_domain "EPOC domain(s)"
label variable intervention_subcategory "Intervention subcategory"
label variable setting "Study setting"
label variable follow_up_str "Follow-up period"
label variable risk_of_bias "Risk of bias"
label variable mortality "Mortality"
label variable hrqol "hrqol"
label variable cost_effectiveness "Cost-effectiveness"
label variable absolute_glycemic_change "Absolute glycemic change"
label variable comparator_description "Comparator description"
label variable provider "Health professional whose behavior changed"

*	Clusters, ICC, and participants variables

label variable total_clusters "Total clusters"
label variable icc "Intraclass coefficient"
label variable total_participants_reported "Total participants enrolled"
label variable part_int_end "Intervention participants, end"
label variable part_con_end "Control participants, end"

*	Mean differences variables

label variable mean_diff_raw "Reported mean difference"
label variable mean_diff_low_ci "Lower 95% CI of mean difference"
label variable mean_diff_upp_ci "Upper 95% CI of mean difference"
label variable mean_diff_p_value "P value of mean difference"

*	Intervention arm variables for a1c

label variable int_mean_change "Within-group change in intervention arm"
label variable int_baseline "Baseline intervention value"
label variable int_baseline_sd "Baseline intervention s.d."
label variable int_baseline_se "Baseline intervention s.e."
label variable int_end "End intervention value"
label variable int_end_sd "End intervention s.d."
label variable int_mean_change_low_ci "Lower 95% CI of within-group change in intervention arm"
label variable int_mean_change_upp_ci "Upper 95% CI of within-group change in intervention arm"
label variable int_sd_raw "Reported s.d. of within-group change in intervention arm"
label variable int_se_raw "Reported s.e. of within-group change in intervention arm"
label variable int_p_value "Intervention change p-value"

*	Control arm variables for a1c

label variable con_mean_change "Within-group change in control arm"
label variable con_baseline "Baseline control value"
label variable con_baseline_sd "Baseline control s.d."
label variable con_baseline_se "Baseline control s.e."
label variable con_end "End control value"
label variable con_end_sd "End control s.d."
label variable con_mean_change_low_ci "Lower 95% CI of within-group change in control arm"
label variable con_mean_change_upp_ci "Upper 95% CI of within-group change in control arm"
label variable con_sd_raw "Reported s.d. of within-group change in control arm"
label variable con_se_raw "Reported s.e. of within-group change in control arm"
label variable con_p_value "Control change p-value"

* Mortality variables
rename con_mortality part_con_death
label variable part_con_death "Control deaths"
rename int_mortality part_int_death
label variable part_int_death "Intervention deaths"

/*******************************************************************************
Part 2: Cleaning data
*******************************************************************************/

//	Sort alphabetically
sort study_title

//	Generating study ID
generate id = _n
label variable id "Study ID"

//	Encoding a few variables that were imported from Excel as string variables

	//	Revising Van Olmen paper's study design to RCT (from separately powered RCT)

		replace study_design = "RCT" if author_year == "Van Olmen 2017"
	
encode study_design, gen(study_design_encoded)
drop study_design

encode risk_of_bias, gen(risk_of_bias_encoded)
drop risk_of_bias

encod epoc_domain, gen(epoc_domain_encoded)
drop epoc_domain

//	Encoding the intervention category manually to allow for post-Excel changes

gen intervention_subcategory_encoded = .
replace intervention_subcategory_encoded = 1 if intervention_subcategory == "Multicomponent clinic-based intervention"
replace intervention_subcategory_encoded = 2 if intervention_subcategory == "Pharmacist-led intervention"
replace intervention_subcategory_encoded = 3 if intervention_subcategory == "Delivery of education, coaching, or support"
replace intervention_subcategory_encoded = 4 if intervention_subcategory == "Case management"
replace intervention_subcategory_encoded = 5 if intervention_subcategory == "Physician training"
replace intervention_subcategory_encoded = 6 if intervention_subcategory == "Multicomponent nurse task shifting"
replace intervention_subcategory_encoded = 7 if intervention_subcategory == "Multicomponent mHealth intervention"
replace intervention_subcategory_encoded = 8 if intervention_subcategory == "Glucose telemonitoring"

label define intervention_subcategory_labels ///
1 "Multicomponent clinic-based intervention"	///
2 "Pharmacist task sharing"	///
3 "Diabetes education or support alone"	///
4 "Case management by nurses"	///
5 "Physician clinical training alone"	///
6 "Multicomponent nurse task sharing"	///
7 "Multicomponent mHealth intervention"	///
8 "Internet-based glucose telemonitoring", modify	//

label values intervention_subcategory_encoded intervention_subcategory_labels
codebook intervention_subcategory_encoded
drop intervention_subcategory
 
//	Adjusting mortality, hrqol, cost effectiveness, a1c change for Table 1 reporting

gen mortality_true = 1
replace mortality_true = 0 if mortality == "NR"

gen hrqol_true = 1
replace hrqol_true = 0 if hrqol == "NR"

gen cost_effectiveness_true = 1
replace cost_effectiveness_true = 0 if cost_effectiveness == "NR"

gen absolute_glycemic_change_true = 1
replace absolute_glycemic_change_true = 0 if absolute_glycemic_change == "NR"

//	Calculating within-group mean differences if not specifically reported

replace int_mean_change = int_end - int_baseline if missing(int_mean_change)

replace con_mean_change = con_end - con_baseline if missing(con_mean_change)

//	Calculating between-group mean differences if not specifically reported

gen mean_diff_calc = mean_diff_raw
replace mean_diff_calc = int_mean_change-con_mean_change if missing(mean_diff_raw)
label variable mean_diff_calc "Calculated mean difference"

/* 	Calculating number of participants who did not experience a death, which is
	necessary for binary meta-analysis in Stata.	*/
	gen part_con_live = .
	replace part_con_live = part_con_end - part_con_death
	
	gen part_int_live = .
	replace part_int_live = part_int_end - part_int_death

/*******************************************************************************
Part 3: Sampling size calculations
*******************************************************************************/	
	
//	Calculate total participants

gen total_participants = .
replace total_participants = part_int_end + part_int_end
label variable total_participants "Total participants in both arms"

//	Calculate design effect = 1 + (M – 1) * ICC
*	See Cochrane link below for more information

gen design_effect = .
replace design_effect = 1 + (total_participants/total_clusters-1)*icc
label variable design_effect "Design effect"
	
//	Calculate effective sample size in intervention and control groups 
* I.e., adjusting for cluster designs for meta-analysis weighting
* https://handbook-5-1.cochrane.org/chapter_16/16_3_4_approximate_analyses_of_cluster_randomized_trials_for_a.htm

codebook study_design_encoded /// Note study design == 1 refers to the cluster RCTs

gen eff_part_int_end = .
replace eff_part_int_end = part_int_end/design_effect if study_design_encoded == 1
replace eff_part_int_end = part_int_end if study_design_encoded == 2
label variable eff_part_int_end "Effective participants in intervention arm"

gen eff_part_con_end = .
replace eff_part_con_end = part_con_end/design_effect if study_design_encoded == 1
replace eff_part_con_end = part_con_end if study_design_encoded == 2
label variable eff_part_con_end "Effective participants in control arm"

//	Stata needs sample size as an integer, so we need to round the effective sample

replace eff_part_int_end = round(eff_part_int_end, 1)
replace eff_part_con_end = round(eff_part_con_end, 1)

/*******************************************************************************
Part 4: Effect size calculations
*******************************************************************************/
	
//	Creating variable of generated intervention and control s.d.

gen int_sd = .
label variable int_sd "Generated s.d. of within-group change in intervention arm"

gen con_sd = .
label variable con_sd "Generated s.d. of within-group change in control arm"

//	Option 1: Using reported s.d. if reported

replace int_sd = int_sd_raw
replace con_sd = con_sd_raw

/*	Option 2 with unreported data: Use within-group confidence intervals, standard 
error, or p-values to derive within-group standard deviations. */
*	See: https://handbook-5-1.cochrane.org/chapter_7/7_7_3_2_obtaining_standard_deviations_from_standard_errors_and.htm */
*	Also see Kirkwood (2003) p. 69 for a reminder of formulas for paired statistics

	*	Option 2a S.d derived from within-group confidence intervals 

	replace int_sd = ((int_mean_change_upp_ci - int_mean_change_low_ci)/3.92)*sqrt(eff_part_int_end) if missing(int_sd)
	replace con_sd = ((con_mean_change_upp_ci - con_mean_change_low_ci)/3.92)*sqrt(eff_part_con_end) if missing(con_sd)

	*	Option 2b: S.d derived from within-group standard errors

	replace int_sd = int_se_raw*sqrt(eff_part_int_end) if missing(int_sd)
	replace con_sd = con_se_raw*sqrt(eff_part_con_end) if missing(con_sd)

	*	Option 2c: S.d derived from within-group p-values
	/*	See a bit more statistical background:
			http://www.stat.columbia.edu/~regina/W1111S09/tutorial5.pdf
			https://handbook-5-1.cochrane.org/chapter_7/7_7_3_3_obtaining_standard_deviations_from_standard_errors.htm	*/
	
	replace int_sd = abs(int_mean_change)/(invttail(eff_part_int_end-1,int_p_value/2)*sqrt(1/eff_part_int_end)) if missing(int_sd)
	replace con_sd = abs(con_mean_change)/(invttail(eff_part_con_end-1,con_p_value/2)*sqrt(1/eff_part_con_end)) if missing(con_sd)
	
/*	Option 3 with unreported data: Use between-group confidence intervals, standard 
error, or p-values to derive within-group standard deviations. */
*	See: https://handbook-5-1.cochrane.org/chapter_7/7_7_3_3_obtaining_standard_deviations_from_standard_errors.htm

	*	*	Option 3a: S.d derived from between-group confidence intervals 

	replace int_sd = ((mean_diff_upp_ci - mean_diff_low_ci)/3.92)/(sqrt(1/eff_part_int_end+1/eff_part_con_end)) if missing(int_sd)

	*	*	Option 3b: S.d derived from between-group p-values  
	
	replace int_sd = (abs(mean_diff_calc)/invttail((eff_part_con_end+eff_part_int_end)-2,mean_diff_p_value/2))/(sqrt(1/eff_part_int_end+1/eff_part_con_end)) if missing(int_sd)

	/*	Note: For all between-group derivations, must assume that s.d. of control 
	group is same as intervention group.
	
	Also: This calculation additionally applies to the Sriram study; this study has
	an inferred p-value of within-group intervention change of 0.01, which allowed
	a within-group s.d. to be inferred for intervention arm. Thus, we are making the assumption
	that, for this study, the s.d. of the control arm equals that of intervention arm. */

	replace con_sd = int_sd if missing(con_sd)	
	
/*	Option 4 with unreported data: Use mean and test statistics of baseline and end
values for each trial arm.	*/
*	See: https://handbook-5-1.cochrane.org/chapter_16/16_1_3_2_imputing_standard_deviations_for_changes_from_baseline.htm

	//	First calculating correlation coefficients available in studies with within-group data
	
	gen corr_int = .
	replace corr_int = (int_baseline_sd^2+int_end_sd^2-int_sd^2)/(2*int_baseline_sd*int_end_sd) ///
		if !missing(int_sd) | !missing(int_mean_change_upp_ci) | !missing(int_se) | !missing(int_p_value)
	
	gen corr_con = .
	replace corr_con = (con_baseline_sd^2+con_end_sd^2-con_sd^2)/(2*con_baseline_sd*con_end_sd) ///
		if !missing(con_sd) | !missing(con_mean_change_upp_ci) | !missing(con_se) | !missing(con_p_value)
	
	summarize corr_int corr_con, detail
	
	/*	Showing that the median correlation within each arm is approximately 0.5
		Since this is also what is often used and recommended in the literature,
		we will use it to impute.	*/
	
	replace int_sd = sqrt(int_baseline_sd^2+int_end_sd^2-(2*0.5*int_baseline_sd*int_end_sd)) if missing(int_sd)
	replace con_sd = sqrt(con_baseline_sd^2+con_end_sd^2-(2*0.5*con_baseline_sd*con_end_sd)) if missing(con_sd)

/*******************************************************************************
Part 5: Combining the interventions in the Anzaldo study
*******************************************************************************/

/*	This bit of code aggregates the 2 intervention arms of Anzaldos study 
compared to the single control arm. */

//	Combined sample size, mean, and s.d.
*	See Cochrane for formula: https://handbook-5-1.cochrane.org/chapter_7/table_7_7_a_formulae_for_combining_groups.htm
/*  Note: I use the below combine function to check calculations using scalars:
https://ideas.repec.org/c/boc/bocode/s457265.html
r(n)                combined group size
r(mean)             combined mean
r(SD)               combined SD

combine n1 mean1 sd1 n2 mean2 sd2		*/

*	Values are entered directly as scalars below
list part_int_end int_mean_change int_sd_raw if author_year == "Anzaldo-Campos 2016, Gilmer 2019" | author_year == "Anzaldo-Campos: Intervention 2"

combine 82 -2.630 3.730 89 -3.020 2.830
 
replace eff_part_int_end = r(n) if author_year == "Anzaldo-Campos 2016, Gilmer 2019"

replace int_mean_change = r(mean) if author_year == "Anzaldo-Campos 2016, Gilmer 2019"

replace int_sd = r(SD) if author_year == "Anzaldo-Campos 2016, Gilmer 2019"

list eff_part_int_end int_mean_change int_sd if author_year == "Anzaldo-Campos 2016, Gilmer 2019" 

//	Drop the superfluous record for Anzaldos

drop if author_year == "Anzaldo-Campos: Intervention 2"

/*******************************************************************************
Part 6: Adjusting the Reutens study
*******************************************************************************/

/*	The Reutens study was a cluster RCT with 10 centers.

LMICs sites:
Beijing - China
Hanoi - Vietnam
Manila - Philippines
Bangkok - Thailand
Surabaya - Indonesia
Johor Bahru - Malaysia.

Non-LMICs sites
Hong Kong
Seoul - South Korea
Taipei - Taiwan
Singapore

Unfortunately, the authors did not respond to email requests as of 2/2. So we will
attempt include 60% of sample into the meta analysis, since 60% of study sites
were in LMICs.
*/

replace eff_part_int_end = eff_part_int_end * 0.6 if author_year == "Reutens 2012"
replace eff_part_con_end = eff_part_con_end * 0.6 if author_year == "Reutens 2012"

//	Need to round the effective sample sizes

replace eff_part_int_end = round(eff_part_int_end, 1) if author_year == "Reutens 2012"
replace eff_part_con_end = round(eff_part_con_end, 1) if author_year == "Reutens 2012"

/*******************************************************************************
Part 7: Double-checking PRISMA diagram
*******************************************************************************/

*	Listing the trials that are missing pre-post HbA1c outcome data
list study_title eff_part_int_end int_mean_change int_sd if risk_of_bias_encoded == 1 | missing(int_sd)

/*******************************************************************************
Part 8: Descriptive statistics
*******************************************************************************/

//	Table 1 info

*	EPOC domain
tabstat total_participants_reported, statistics( sum count) by(epoc_domain_encoded)

*	Intervention type
tabulate intervention_subcategory_encoded, missing

*	Study design
tabulate study_design_encoded, missing

*	World Bank region
tabulate wb_region, missing

*	World Bank income groups
tabulate wb_income_group, missing

*	Study setting
tabulate setting, missing

*	Outcomes reported	
tabstat mortality_true hrqol_true cost_effectiveness_true absolute_glycemic_change_true, statistics( sum ) by(intervention_subcategory_encoded)

*	Study duration and follow-up

//	The duration and follow-up variable needs to be converted to a number and destringed
	replace intervention_duration = "12 months" if intervention_duration == "1 year"
	replace intervention_duration = "24 months" if intervention_duration == "2 years"
	ereplace intervention_duration = sieve(intervention_duration), keep(numeric)
	destring intervention_duration, replace
	
	replace follow_up_str = "12 months" if follow_up_str == "1 year"
	replace follow_up_str = "24 months" if follow_up_str == "2 years"
	ereplace follow_up_str = sieve(follow_up_str), keep(numeric)
	destring follow_up_str, replace
	replace follow_up_str = 12.5 if author_year == "Tutino 2017"
	
	tabstat intervention_duration, statistics( min p25 median p75 max )
	tabstat follow_up_str, statistics( min p25 median p75 max )
	
*	Risk of bias

tab risk_of_bias_encoded
codebook risk_of_bias_encoded
list study_title if risk_of_bias_encoded == 2

*	Making excel table of study information, Table S1
export excel study_title other_reports epoc_domain intervention_subcategory   country setting_string intervention_description intervention_duration follow_up_str provider mortality hrqol cost_effectiveness absolute_glycemic_change using "/Users/davidflood/Desktop/Table S1.xlsx", firstrow(variables) replace

/*******************************************************************************
Part 9: GRADE assessment information
*******************************************************************************/

tabstat total_participants_reported, statistics( sum ) by(intervention_subcategory_encoded)

/*******************************************************************************
Part 10: Conducting the meta-analysis
*******************************************************************************/

* Review studies included in meta analysis  ------------------------------------

tab intervention_subcategory_encoded epoc_domain_encoded
tab study_title epoc_domain_encoded
tab study_title intervention_subcategory_encoded

format int_baseline con_baseline int_mean_change int_sd con_mean_change con_sd mean_diff_raw mean_diff_calc %3.1f

list study_title int_baseline eff_part_int_end eff_part_con_end con_baseline int_mean_change int_sd con_mean_change con_sd mean_diff_raw mean_diff_calc, string(5)

list study_title epoc_domain_encoded intervention_subcategory_encoded, string(15)
tab intervention_subcategory_encoded epoc_domain_encoded

* Forest plot for mortality ----------------------------------------------------

meta esize part_int_death part_int_live part_con_death part_con_live, random(dlaird) studylabel(study_title) 

meta summarize, cformat(%9.2f) subgroup(intervention_subcategory_encoded) eform

meta forest, subgroup(intervention_subcategory_encoded) name(mortality_forest, replace) cibind(parentheses) markeropts(mcolor(gs10))  insidemarker(mcolor(black) msize(tiny) msymbol(diamond))  nonotes ciopts(lcolor(black)) nullrefline(favorsright("Favors comparator", margin(medium)) favorsleft("Favors intervention`'", margin(large)) lcolor(black))   esrefline(lcolor(black) lpattern(vshortdash)) noohetstats noohomtest noosigtest noghetstats nogwhomtests  nogbhomtests  nogmarkers nooverall eform

	* Formatting via graph edit
{
	gr_edit .plotregion1.column9.group_items[1].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	// group_items[1] edits

	gr_edit .plotregion1.column10.group_items[1].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	// group_items[1] edits

	gr_edit .plotregion1.column11.group_items[1].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	// group_items[1] edits

	gr_edit .plotregion1.column9.group_items[2].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	// group_items[2] edits

	gr_edit .plotregion1.column10.group_items[2].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	// group_items[2] edits

	gr_edit .plotregion1.column11.group_items[2].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	// group_items[2] edits

	gr_edit .plotregion1.column9.group_items[3].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	// group_items[3] edits

	gr_edit .plotregion1.column10.group_items[3].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	// group_items[3] edits

	gr_edit .plotregion1.column11.group_items[3].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	// group_items[3] edits

	gr_edit .plotregion1.column9.group_items[4].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	// group_items[4] edits

	gr_edit .plotregion1.column10.group_items[4].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	// group_items[4] edits

	gr_edit .plotregion1.column11.group_items[4].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	// group_items[4] edits

	gr_edit .plotregion1.column9.group_items[5].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	// group_items[5] edits

	gr_edit .plotregion1.column10.group_items[5].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	// group_items[5] edits

	gr_edit .plotregion1.column11.group_items[5].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	// group_items[5] edits

	gr_edit .plotregion1.column9.group_items[6].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	// group_items[6] edits

	gr_edit .plotregion1.column10.group_items[6].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	// group_items[6] edits

	gr_edit .plotregion1.column11.group_items[6].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	// group_items[6] edits

	// <end>
}

graph save "Graphs/mortality_forest.gph", replace
graph export "Graphs/mortality_forest.eps", as(eps) name("mortality_forest") preview(on) replace

* Primary meta-analysis for HbA1c ---------------------------------------------

//	Set up meta analysis function and generate forest plot for full analysis
* 	Note syntax: meta esize nt meant sdt nc meanc sdc

codebook epoc_domain_encoded

meta esize eff_part_int_end int_mean_change int_sd eff_part_con_end con_mean_change con_sd, esize(mdiff) random(dlaird) studylabel(study_title)

meta summarize, cformat(%9.2f) predinterval    // to generate prediction interval
meta summarize, cformat(%9.2f)subgroup(intervention_subcategory_encoded)   

meta forest, columnopts(_weight, format(%9.1f)) columnopts(_mean1 _sd1 _mean2 _sd2, format(%9.2f))	///
	subgroup(intervention_subcategory_encoded) 						///
	name(main_forest, replace) 										///
	markeropts(mcolor(gs10)) 										///
	gmarkeropts(mcolor(black)  										///
		mfcolor(none))  											///
	insidemarker(mcolor(black) 										///
		msize(tiny) msymbol(diamond)) 								///
	nonotes 														///
	ciopts(lcolor(black)) 											///
	nullrefline(favorsright("Favors comparator", margin(medium)) favorsleft("Favors intervention", margin(large)) lcolor(black)) ///
	omarkeropts(mfcolor(black) mcolor(black)) 						///
	esrefline(lcolor(black) lpattern(vshortdash))  					///
	nogwhomtests 													///
		nogbhomtests 												//

{
	// Plot grec edits
	// hom_stats[1] edits
	gr_edit .plotregion1.column1.hom_stats[1].text = {}
	gr_edit .plotregion1.column1.hom_stats[1].text.Arrpush Heterogeneity: {&tau}{sup:2} = 0.13, I{sup:2} = 87.8% [84.0, 90.6], p<0.001
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	
	// het_stats[9] edits
	gr_edit .plotregion1.column1.het_stats[9].text = {}
	gr_edit .plotregion1.column1.het_stats[9].text.Arrpush Prediction interval: [-1.19, 0.28]
	gr_edit .set_spositions
	gr_edit .set_graphwidth

	// het_stats[1] edits
	gr_edit .plotregion1.column1.het_stats[1].text = {}
	gr_edit .plotregion1.column1.het_stats[1].text.Arrpush Heterogeneity: {&tau}{sup:2} = 0.06, I{sup:2} = 60.0% [8.2%, 82.6%], p=0.020
	gr_edit .set_spositions
	gr_edit .set_graphwidth

	// het_stats[2] edits
	gr_edit .plotregion1.column1.het_stats[2].text = {}
	gr_edit .plotregion1.column1.het_stats[2].text.Arrpush Heterogeneity: {&tau}{sup:2} = 0.29, I{sup:2} = 91.0% [86.5, 94.0], p<0.001
	gr_edit .set_spositions
	gr_edit .set_graphwidth

	// het_stats[3] edits
	gr_edit .plotregion1.column1.het_stats[3].text = {}
	gr_edit .plotregion1.column1.het_stats[3].text.Arrpush Heterogeneity: {&tau}{sup:2} = 0.06, I{sup:2} = 64.1% [18.8, 84.1], p=0.010
	gr_edit .set_spositions
	gr_edit .set_graphwidth

	// het_stats[4] edits
	gr_edit .plotregion1.column1.het_stats[4].text = {}
	gr_edit .plotregion1.column1.het_stats[4].text.Arrpush Heterogeneity: {&tau}{sup:2} = 0.03, I{sup:2} = 46.8%, p=0.170
	gr_edit .set_spositions
	gr_edit .set_graphwidth

	// het_stats[5] edits
	gr_edit .plotregion1.column1.het_stats[5].text = {}
	gr_edit .plotregion1.column1.het_stats[5].text.Arrpush Heterogeneity: {&tau}{sup:2} = 0.02, I{sup:2} = 41.4%, p=0.191
	gr_edit .set_spositions
	gr_edit .set_graphwidth

	// het_stats[6] edits
	gr_edit .plotregion1.column1.het_stats[6].text = {}
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	
	// het_stats[7] edits
	gr_edit .plotregion1.column1.het_stats[7].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	
		// het_stats[8] edits
	gr_edit .plotregion1.column1.het_stats[8].text = {}
	gr_edit .plotregion1.column1.het_stats[8].text.Arrpush Heterogeneity: {&tau}{sup:2} = 0.15, I{sup:2} = 93.8%, p<0.001
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	
	// group marker[6] edits
	gr_edit .plotregion1.column9.group[6].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth

	// group_items[6] edits
	gr_edit .plotregion1.column11.group_items[6].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth

	// group_items[6] edits
	gr_edit .plotregion1.column12.group_items[6].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth

	// group_items[6] edits
	gr_edit .plotregion1.column13.group_items[6].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth

	// group marker[7] edits
	gr_edit .plotregion1.column9.group[7].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth

	// group_items[7] edits
	gr_edit .plotregion1.column11.group_items[7].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth

	// group_items[7] edits
	gr_edit .plotregion1.column12.group_items[7].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth

	// group_items[7] edits
	gr_edit .plotregion1.column13.group_items[7].draw_view.setstyle, style(no)
	gr_edit .set_spositions
	gr_edit .set_graphwidth

	// line[1] edits
	gr_edit .plotregion1.column9.AddLine added_lines editor -1.258580963413391 90.52975321221979 .178556789420852 90.64594329933594
	gr_edit .plotregion1.column9.added_lines_new = 1
	gr_edit .plotregion1.column9.added_lines_rec = 1
	gr_edit .plotregion1.column9.added_lines[1].style.editstyle  linestyle( width( sztype(relative) val(.2) allow_pct(1)) color(black) pattern(solid) align(inside)) headstyle( symbol(circle) linestyle( width( sztype(relative) val(.2) allow_pct(1)) color(black) pattern(solid) align(inside)) fillcolor(black) size( sztype(relative) val(1.52778) allow_pct(1)) angle(stdarrow) symangle(zero) backsymbol(none) backline( width( sztype(relative) val(.2) allow_pct(1)) color(black) pattern(solid) align(inside)) backcolor(black) backsize( sztype(relative) val(0) allow_pct(1)) backangle(stdarrow) backsymangle(zero)) headpos(neither) editcopy
	gr_edit .plotregion1.column9.added_lines[1].drag_point = (0)
	gr_edit .plotregion1.column9.added_lines[1].DragBy .3485702613484613 .1596819725371406
	gr_edit .plotregion1.column9.added_lines[1].drag_point = (0)
	gr_edit .plotregion1.column9.added_lines[1].DragBy -.1161900871161523 0
	gr_edit .plotregion1.column9.added_lines[1].style.editstyle linestyle(width(thick)) editcopy
	gr_edit .plotregion1.column9.added_lines[1].style.editstyle linestyle(pattern(solid)) editcopy
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	gr_edit .plotregion1.column9.added_lines[1].drag_point = (0)
	gr_edit .plotregion1.column9.added_lines[1].DragBy -.2323801742323054 -.0798409862685752
	gr_edit .plotregion1.column9.added_lines[1].drag_point = (0)
	gr_edit .plotregion1.column9.added_lines[1].DragBy .5809504355807639 1.516978739102795
	gr_edit .plotregion1.column9.added_lines[1].drag_point = (0)
	gr_edit .plotregion1.column9.added_lines[1].DragBy -.5809504355807639 -1.676660711639935
	gr_edit .plotregion1.column9.added_lines[1].drag_point = (0)
	gr_edit .plotregion1.column9.added_lines[1].DragBy .3485702613484613 .0798409862685654
	gr_edit .plotregion1.column9.added_lines[1].drag_point = (0)
	gr_edit .plotregion1.column9.added_lines[1].DragBy -.3485702613484578 0
	gr_edit .plotregion1.column9.added_lines[1].drag_point = (0)
	gr_edit .plotregion1.column9.added_lines[1].DragBy .1161900871161523 0
	
		
	gr_edit .plotregion1.column2.AddTextBox added_text editor 90.52973880198901 .3309286875055787
	gr_edit .plotregion1.column2.added_text_new = 1
	gr_edit .plotregion1.column2.added_text_rec = 1
	gr_edit .plotregion1.column2.added_text[1].style.editstyle  angle(default) size( sztype(relative) val(3.4722) allow_pct(1)) color(black) horizontal(left) vertical(middle) margin( gleft( sztype(relative) val(0) allow_pct(1)) gright( sztype(relative) val(0) allow_pct(1)) gtop( sztype(relative) val(0) allow_pct(1)) gbottom( sztype(relative) val(0) allow_pct(1))) linegap( sztype(relative) val(0) allow_pct(1)) drawbox(no) boxmargin( gleft( sztype(relative) val(0) allow_pct(1)) gright( sztype(relative) val(0) allow_pct(1)) gtop( sztype(relative) val(0) allow_pct(1)) gbottom( sztype(relative) val(0) allow_pct(1))) fillcolor(bluishgray) linestyle( width( sztype(relative) val(.2) allow_pct(1)) color(black) pattern(solid) align(inside)) box_alignment(east) editcopy
	gr_edit .plotregion1.column2.added_text[1].text = {}
	gr_edit .plotregion1.column2.added_text[1].text.Arrpush {bf:5,240}
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	// editor text[1] edits

	gr_edit .plotregion1.column5.AddTextBox added_text editor 90.64592915949537 .1773967100258255
	gr_edit .plotregion1.column5.added_text_new = 1
	gr_edit .plotregion1.column5.added_text_rec = 1
	gr_edit .plotregion1.column5.added_text[1].style.editstyle  angle(default) size( sztype(relative) val(3.4722) allow_pct(1)) color(black) horizontal(left) vertical(middle) margin( gleft( sztype(relative) val(0) allow_pct(1)) gright( sztype(relative) val(0) allow_pct(1)) gtop( sztype(relative) val(0) allow_pct(1)) gbottom( sztype(relative) val(0) allow_pct(1))) linegap( sztype(relative) val(0) allow_pct(1)) drawbox(no) boxmargin( gleft( sztype(relative) val(0) allow_pct(1)) gright( sztype(relative) val(0) allow_pct(1)) gtop( sztype(relative) val(0) allow_pct(1)) gbottom( sztype(relative) val(0) allow_pct(1))) fillcolor(bluishgray) linestyle( width( sztype(relative) val(.2) allow_pct(1)) color(black) pattern(solid) align(inside)) box_alignment(east) editcopy
	gr_edit .plotregion1.column5.added_text[1].text = {}
	gr_edit .plotregion1.column5.added_text[1].text.Arrpush {bf:4,820}
	gr_edit .set_spositions
	gr_edit .set_graphwidth
	// editor text[1] edits

	gr_edit .plotregion1.column5.added_text[1].DragBy -.1161903575063788 -.1180983578263724
	// editor text[1] reposition

	gr_edit .plotregion1.column2.added_text[1].DragBy .2323807150127569 -.3149289542036615
	// editor text[1] reposition

	gr_edit .plotregion1.column2.added_text[1].DragBy -1.045713217557409 .2952458945659329
	// editor text[1] reposition

	gr_edit .plotregion1.column2.added_text[1].DragBy .5809517875318922 -.2558797752904758
	// editor text[1] reposition

	gr_edit .plotregion1.column2.added_text[1].DragBy .3485710725191387 -.0196830596377291
	// editor text[1] reposition

	gr_edit .plotregion1.column2.added_text[1].DragBy -.1161903575063788 0
	// editor text[1] reposition

}

graph save "Graphs/main_forest.gph", replace		
graph export "Graphs/main_forest.eps", as(eps) name("main_forest") preview(on) replace

//	Subtotals for forest, to be added manually
tab study_title if !missing(_meta_es)

tabstat eff_part_int_end eff_part_con_end if !missing(_meta_es) , statistics( sum ) by(intervention_subcategory_encoded)


******************** SENSITIVITY AND HETEROGENEITY *****************************

/*	Note: The best way to edit the forest plot graphic is as follows:
1.	Open gph file in Stata on Windows (can use a virtual network at umich)
2. 	Save as enchanged meta file (EMF)
3.	Open the EMF file in Powerpoint/Office on Windows and ungroup the graphic
4.	Save this Powerpoint/Office file and now can open in Mac and revise as needed	*/

meta esize eff_part_int_end int_mean_change int_sd eff_part_con_end con_mean_change con_sd, esize(mdiff) random(dlaird) studylabel(study_title)

* Trim and fill ----------------------------------------------------------------

preserve

drop if inlist(study_title,"Chao","Khetan","Zhong","Gillani")
	// For some reason you need to drop the studies not included in the analysis to get trim and fill to work

meta esize eff_part_int_end int_mean_change int_sd eff_part_con_end con_mean_change con_sd, esize(mdiff) random(dlaird) studylabel(study_title)

meta trimfill, funnel(scheme(s1mono) title(Trim-and-fill funnel plot) xtitle(Between-arm mean HbA1c change (%)) ytitle("Standard error")) itermethod() poolmethod() log

graph save "Graphs/trim_fill.gph", replace
graph export "Graphs/trim_fill.eps", as(eps) name("Graph") preview(on) replace

restore

*	Funnel plot and Egger overall
meta funnelplot, name(funnel_main, replace) scheme(s1mono) note("Egger p = <0.001") legend(off) title("Overall funnel plot (all studies)") xtitle("Between-arm mean HbA1c change (%)") ytitle("Standard error")  random(dlaird)

graph save "Graphs/funnel_main.gph", replace
graph export "Graphs/funnel_main.eps", as(eps) name("funnel_main") preview(on) replace

meta bias, egger detail

*	Funnel plot and Egger by intervention type  --------------------------------

//	Calculate Egger p for subcategories
	foreach x of numlist 1 2 3 4 5 8 {
	preserve
	drop if intervention_subcategory_encoded != `x'
	tab intervention_subcategory_encoded
	meta bias, egger detail
	restore
	}

meta funnelplot, name(funnel_subcategory, replace) by(intervention_subcategory_encoded, note("") legend(off) title("Funnel plot by intervention type (all studies)", size(*0.9))) scheme(s1mono) xtitle("Between-arm mean HbA1c change (%)") ytitle("Standard error")  subtitle(,size(*0.7)) random(dlaird) xscale(range(-3.2 2.2)) xlabel (-3(1)2) yscale(range(-0.1 1.1)) ylabel (0(0.5)1)

	//	Code to edit graph and include Egger p values
		gr_edit .plotregion1.xaxis1[1].draw_view.setstyle, style(yes)
		// xaxis1[1] edits

		gr_edit .plotregion1.xaxis1[2].draw_view.setstyle, style(yes)
		// xaxis1[2] edits

		gr_edit .plotregion1.xaxis1[3].draw_view.setstyle, style(yes)
		// xaxis1[3] edits

		gr_edit .plotregion1.xaxis1[4].draw_view.setstyle, style(yes)
		// xaxis1[4] edits

		gr_edit .plotregion1.xaxis1[5].draw_view.setstyle, style(yes)
		// xaxis1[5] edits

		gr_edit .plotregion1.legend[1].draw_view.setstyle, style(yes)
		// legend[1] edits

		gr_edit .plotregion1.legend[1].draw_view.setstyle, style(no)
		// legend[1] edits

		gr_edit .plotregion1.note[1].text = {}
		gr_edit .plotregion1.note[1].text.Arrpush Egger p = 0.15
		// note[1] edits

		gr_edit .plotregion1.note[2].text = {}
		gr_edit .plotregion1.note[2].text.Arrpush Egger p = 0.13
		// note[2] edits

		gr_edit .plotregion1.note[3].text = {}
		gr_edit .plotregion1.note[3].text.Arrpush Egger p = 0.12
		// note[3] edits
		
graph save "Graphs/funnel_subcategory.gph", replace	
graph export "Graphs/funnel_subcategory.eps", as(eps) name("funnel_subcategory") preview(on) replace
	
/*	Note: This code runs the meta bias function only on the 5 intervention subcategories
in which there are >1 study to pool.	*/

codebook intervention_subcategory_encoded

* Leave-one-out meta-analysis --------------------------------------------------

/* 	This is the leave-one-out method, known generally as a set of methods of "influence analysis"

This is an excellent summary in R:
https://bookdown.org/MathiasHarrer/Doing_Meta_Analysis_in_R/influenceanalyses.html

In Stata the admetan package can do this analysis but not the new native meta package.

You can confirm the underlying estimates using the following code, which are equivalent
to those used in the native meta esize package:

metan eff_part_int_end int_mean_change int_sd eff_part_con_end con_mean_change con_sd, randomi label(namevar=study_title) nostandard

It is useful to sort the leave-one-out results to facilitate interpretation. 	*/

admetan eff_part_int_end int_mean_change int_sd eff_part_con_end con_mean_change con_sd, model(dl) label(namevar=study_title) nostandard influence nograph  // need to use nograph function to permit sortby in the following line

admetan eff_part_int_end int_mean_change int_sd eff_part_con_end con_mean_change con_sd, model(dl) label(namevar=study_title) nostandard influence sortby(_ES) forestplot( xlabel(-1(0.5)0) name("leave_one_out",replace)) effect(Mean diff.)

	* Formatting via graph edit
{	
	// Graph edits
	gr_edit .style.editstyle margin(zero) editcopy
	gr_edit .style.editstyle boxstyle(shadestyle(color(white))) editcopy
	gr_edit .style.editstyle boxstyle(linestyle(color(white))) editcopy
	
	// plot9 edits
	gr_edit .plotregion1.plot9.draw_view.setstyle, style(no)
	
	// Graph size
	gr_edit .style.editstyle declared_ysize(8) editcopy
	gr_edit .style.editstyle declared_xsize(8) editcopy
}

graph save "Graphs/leave_one_out.gph", replace	
graph export "Graphs/leave_one_out.eps", as(eps) name(leave_one_out) preview(on) replace
	
/* This code then removes the Javaid study to assess what the I2 is when this
study is excluded.	*/
	
meta esize eff_part_int_end int_mean_change int_sd eff_part_con_end con_mean_change con_sd if study_title !="Javaid", esize(mdiff) random(dlaird) studylabel(study_title)

meta summarize, cformat(%9.2f) subgroup(intervention_subcategory_encoded)

* Generating forest plot with high risk of bias studies excluded  --------------

//	Resetring meta esize

meta esize eff_part_int_end int_mean_change int_sd eff_part_con_end con_mean_change con_sd if risk_of_bias_encoded !=1, esize(mdiff) random(dlaird) studylabel(study_title)

meta forest if risk_of_bias_encoded !=1, columnopts(_mean1 _sd1 _mean2 _sd2, format(%9.2f)) subgroup(intervention_subcategory_encoded) name(sensitivity_forest, replace) cibind(parentheses) markeropts(mcolor(gs10)) gmarkeropts(mcolor(black) mfcolor(none)) insidemarker(mcolor(black) msize(tiny) msymbol(diamond))  nonotes ciopts(lcolor(black)) nullrefline(favorsright("Favors comparator", margin(medium)) favorsleft("Favors intervention", margin(large)) lcolor(black)) noohomtest nogwhomtests nogbhomtests  omarkeropts(mfcolor(black) mcolor(black))  esrefline(lcolor(black) lpattern(vshortdash))

graph save "Graphs/sensitivity_forest.gph", replace
graph export "Graphs/sensitivity_forest.eps", as(eps) name("sensitivity_forest") preview(on) replace

*	Funnel plot and Egger overall, restricting high risk of bias studies

meta funnelplot if risk_of_bias_encoded !=1, name(funnel_main_rob, replace)  scheme(s1mono) note("Egger p = 0.002") legend(off) title("Overall funnel plot (excluding studies at high risk of bias)") xtitle("Between-arm mean HbA1c change (%)") ytitle("Standard error") random(dlaird)
	
graph save "Graphs/funnel_main_rob.gph", replace
graph export "Graphs/funnel_main_rob.eps", as(eps) name("funnel_main_rob") preview(on) replace

meta bias if risk_of_bias_encoded !=1, egger detail

*	Funnel plot and Egger by intervention type, restricting high risk of bias studies
	
	//	Calculate Egger p for subcategories
	foreach x of numlist 1 2 3 4 {
	preserve
	drop if intervention_subcategory_encoded != `x'
	meta bias if risk_of_bias_encoded !=1, egger detail
	restore
}

meta funnelplot if risk_of_bias_encoded !=1, name(funnel_subcategory_rob, replace) by(intervention_subcategory_encoded, note("") legend(off) title("Funnel plot by intervention type (excluding studies at high risk of bias)", size(*0.8))) scheme(s1mono) xtitle("Between-arm mean HbA1c change (%)") ytitle("Standard error")  subtitle(,size(*0.7))  note("Egger p = ") random(dlaird) xscale(range(-3.2 2.2)) xlabel (-3(1)2) yscale(range(-0.1 1.1)) ylabel (0(0.5)1)

	//	Code to edit graph and include Egger p values
{
	gr_edit .plotregion1.yaxis1[8].draw_view.setstyle, style(yes)
	// yaxis1[8] edits

	gr_edit .plotregion1.xaxis1[2].draw_view.setstyle, style(yes)
	// xaxis1[2] edits

	gr_edit .plotregion1.xaxis1[1].draw_view.setstyle, style(yes)
	// xaxis1[1] edits

	gr_edit .plotregion1.xaxis1[3].draw_view.setstyle, style(yes)
	// xaxis1[3] edits

	gr_edit .plotregion1.xaxis1[4].draw_view.setstyle, style(yes)
	// xaxis1[4] edits

	gr_edit .plotregion1.xaxis1[5].draw_view.setstyle, style(yes)
	// xaxis1[5] edits

	gr_edit .plotregion1.note[4].draw_view.setstyle, style(no)
	// note[4] edits

	gr_edit .plotregion1.note[5].draw_view.setstyle, style(no)
	// note[5] edits

	gr_edit .plotregion1.note[6].draw_view.setstyle, style(no)
	// note[6] edits

	gr_edit .plotregion1.note[7].draw_view.setstyle, style(no)
	// note[7] edits

	gr_edit .plotregion1.note[8].draw_view.setstyle, style(no)
	// note[8] edits

	gr_edit .plotregion1.note[1].text = {}
	gr_edit .plotregion1.note[1].text.Arrpush Egger p = 0.15
	// note[1] edits

	gr_edit  .plotregion1.note[2].text = {}
	gr_edit  .plotregion1.note[2].text.Arrpush Egger p = 0.10
	// note[2] edits

	gr_edit .plotregion1.note[3].text = {}
	gr_edit .plotregion1.note[3].text.Arrpush Egger p = 0.36
	// note[3] edits
}


graph save "Graphs/funnel_subcategory_rob.gph", replace
graph export "Graphs/funnel_subcategory_rob.eps", as(eps) name("funnel_subcategory_rob") preview(on) replace
	
tab intervention_subcategory_encoded if risk_of_bias_encoded !=1
codebook intervention_subcategory_encoded



/*******************************************************************************
Generating the I2 confidence intervals
*******************************************************************************/

* Making an I2 program ---------------------------------------------------------

/* See p. 124 of the following reference for the formulas:
	
	Borenstein M, Hedges LV, Higgins JP, Rothstein HR. Introduction to meta-analysis.
	John Wiley & Sons; 2011.	*/
	
program drop _all
program define i2_confidence_interval 

	if Q > df + 1 {
		scalar B = 0.5*((ln(Q)-ln(df))/(sqrt(2*Q)-sqrt(2*df-1)))
		}
	
	else {
		scalar B = sqrt(1/(2*(df-1)*(1-(1/(3*(df-1)^2)))))
	}
	
	di "B = " B

	scalar L = exp(0.5*ln(Q/df)-1.96*B)
	di "L = " L
	
	scalar U = exp(0.5*ln(Q/df)+1.96*B)
	di "U = " U

	scalar I2 = ((Q-df)/Q)*100
	di "I2 = " I2
	
	scalar lower_ci_i2 = ((L^2-1)/L^2)*100
	di "lower bound of 95% CI = " lower_ci_i2

	scalar upper_ci_i2 = ((U^2-1)/U^2)*100
	di "upper bound of 95% CI = " upper_ci_i2
end	

* Generating the confidence intervals for I2 in the main analysis --------------

meta esize eff_part_int_end int_mean_change int_sd eff_part_con_end con_mean_change con_sd, esize(mdiff) random(dlaird) studylabel(study_title)

meta summarize, cformat(%9.2f) subgroup(intervention_subcategory_encoded)



	* Overall
		scalar Q = 278.21
		scalar df = 34	
		i2_confidence_interval

	* Multicomponent clinic-based intervention
		scalar Q = 15
		scalar df = 6
		i2_confidence_interval

	* Pharmacist-led intervention
		scalar Q = 133.79
		scalar df = 12	
		i2_confidence_interval

	* Delivery of education, coaching, or support
		scalar Q = 16.69 
		scalar df = 6	
		i2_confidence_interval
		
	* Case management
		scalar Q = 1.88
		scalar df = 1	
		i2_confidence_interval

	* Physician training
		scalar Q = 1.71
		scalar df = 1
		i2_confidence_interval

	* Glucose telemonitoring
		scalar Q = 16.15
		scalar df = 1
		i2_confidence_interval

* For the overall meta-analysis excluding studies with high risk of bias -------

meta esize eff_part_int_end int_mean_change int_sd eff_part_con_end con_mean_change con_sd if risk_of_bias_encoded !=1, esize(mdiff) random(dlaird) studylabel(study_title)

meta summarize, cformat(%9.2f) subgroup(intervention_subcategory_encoded)

	* Overall
		scalar Q = 70.83
		scalar df = 20
		i2_confidence_interval

	* Multicomponent clinic-based intervention
		scalar Q = 15
		scalar df = 6
		i2_confidence_interval

	* Pharmacist-led intervention
		scalar Q = 28.79
		scalar df = 3
		i2_confidence_interval

	* Delivery of education, coaching, or support
		scalar Q = 1.5
		scalar df = 3
		i2_confidence_interval
	
/*******************************************************************************
Other random analyses
*******************************************************************************/	

// Proportion of studies with HbA1c as primary outcome
/* Note: Following 3 studies contain HbA1c but for powering:
Gillani 2016 (telemonitoring + control)
Goruntla 2019
Paz-Pacheco 2017		*/

list study_title primary_outcome

gen hba1c_outcome = strpos(primary_outcome, "HbA1c") > 0
tab hba1c_outcome

// For S15 Appendix: Summary of findings table
tabstat total_participants_reported, statistics( sum count) by(intervention_subcategory_encoded)
